── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract() masks raster::extract()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
✖ dplyr::select() masks raster::select()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Attaching package: 'scales'
The following object is masked from 'package:purrr':
discard
The following object is masked from 'package:readr':
col_factor
library(sf)
Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(gt)library(patchwork)
Attaching package: 'patchwork'
The following object is masked from 'package:raster':
area
library(tidyquant)
Loading required package: PerformanceAnalytics
Loading required package: xts
Loading required package: zoo
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
######################### Warning from 'xts' package ##########################
# #
# The dplyr lag() function breaks how base R's lag() function is supposed to #
# work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
# source() into this session won't work correctly. #
# #
# Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
# conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
# dplyr from breaking base R's lag() function. #
# #
# Code in packages is not affected. It's protected by R's namespace mechanism #
# Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
# #
###############################################################################
Attaching package: 'xts'
The following objects are masked from 'package:dplyr':
first, last
Attaching package: 'PerformanceAnalytics'
The following object is masked from 'package:graphics':
legend
Loading required package: quantmod
Loading required package: TTR
Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo
conflicted::conflict_prefer("select", "dplyr")
[conflicted] Will prefer dplyr::select over any other package.
conflicted::conflict_prefer("filter", "dplyr")
[conflicted] Will prefer dplyr::filter over any other package.
# Plot themetheme_set(theme_gmri(rect =element_rect(fill ="white", color =NA), legend.position ="bottom"))# Degree symboldeg_c <-"\u00b0C"# vectors for factor levelsarea_levels <-c("GoM", "GB", "SNE", "MAB")area_levels_long <-c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")
Rows: 2976 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): survey_area, month
dbl (2): bot_temp, year
date (1): date
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
`summarise()` has grouped output by 'survey_area'. You can override using the
`.groups` argument.
Joining with `by = join_by(survey_area, month)`
# Load polygons for aggregate regions based on trawl survey strataregions <-read_sf(here::here("Data/raw", "nmfs_trawl_regions_collection.geojson"))# Prepare the regions as single polygons:gb <-filter(regions, area =="Georges Bank") %>%st_union() %>%st_as_sf()gom <-filter(regions, area =="Gulf of Maine") %>%st_union() %>%st_as_sf()mab <-filter(regions, area =="Mid-Atlantic Bight") %>%st_union() %>%st_as_sf()sne <-filter(regions, area =="Southern New England") %>%st_union() %>%st_as_sf()
# month weighting so we don't have to go back to daily when averagingseas_month_wts <-data.frame("month"=c("03", "04", "05", "09", "10", "11"),"n_days"=days_in_month(as.Date(c("2000-03-01","2000-04-01","2000-05-01","2000-09-01","2000-10-01","2000-11-01"))))# Calculate average for the seasonseasonal_temps <- regional_clim %>%left_join(seas_month_wts) %>%mutate(season =ifelse( month %in%str_pad(c(3:5), side ="left", pad ="0", width =2), "Spring", NA),season =ifelse( month %in%str_pad(c(9:11), side ="left", pad ="0", width =2), "Fall", season),survey_area =factor(survey_area, levels = area_levels_long)) %>%drop_na(season) %>%group_by(survey_area, year, season) %>%summarise(sst =weighted.mean(sst, w = n_days),clim_sst =weighted.mean(clim_sst, w = n_days),sst_anom =weighted.mean(sst_anom, w = n_days),bot_temp =weighted.mean(bot_temp, w = n_days),clim_bt =weighted.mean(clim_bt, w = n_days),bt_anom =weighted.mean(bt_anom, w = n_days),.groups ="drop" )
Joining with `by = join_by(month)`
# Calculate annual averagesannual_temps <- regional_clim %>%mutate(n_days =days_in_month(date)) %>%group_by(survey_area, year) %>%summarise(sst =weighted.mean(sst, w = n_days),clim_sst =weighted.mean(clim_sst, w = n_days),sst_anom =weighted.mean(sst_anom, w = n_days),bot_temp =weighted.mean(bot_temp, w = n_days),clim_bt =weighted.mean(clim_bt, w = n_days),bt_anom =weighted.mean(bt_anom, w = n_days),.groups ="drop" )
Timeseries of Monthly Surface and Bottom Temperatures
The following plot shows what the monthly average temperatures have been based on these two datasets.
ggplot(regional_clim) +geom_line(aes(date, sst, color ="SST")) +geom_line(aes(date, bot_temp, color ="BT")) +scale_color_gmri(reverse = T) +facet_wrap(~survey_area, ncol =1) +labs(y ="Temperature", color ="Color")
Warning: Removed 272 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 45 rows containing missing values or values outside the scale range
(`geom_line()`).
Timeseries of 1991-2020 Anomalies
Both surface and bottom span this 30-year period so we can view them as anomalies based on the same reference period.
.{panel-tabset}
Annual
ggplot(regional_clim) +geom_line(aes(date, sst_anom, color ="Surface Temperature Anomalies")) +geom_line(aes(date, bt_anom, color ="Bottom Temperature Anomalies")) +scale_color_gmri(reverse = T) +facet_wrap(~survey_area, ncol =1) +labs(y ="Temperature Anomaly", color ="Color")
Warning: Removed 272 rows containing missing values or values outside the scale range
(`geom_line()`).
Warning: Removed 45 rows containing missing values or values outside the scale range
(`geom_line()`).
Seasonal
# Make these seasonal as if we were matching it to seasonal trawl data# Calculate average for the seasonseasonal_temps %>%rename(`Bottom Temperature`= bt_anom, `Surface Temperature`= sst_anom) %>%pivot_longer(names_to ="var", values_to ="vals", cols =ends_with("Temperature")) %>%ggplot(aes(year, vals, color = survey_area)) +geom_point(alpha =0.35, size =0.5) +geom_ma(aes(linetype ="5-Year Moving Average"), n =5, ma_fun = SMA) +facet_grid(fct_rev(var) ~fct_rev(season)) +scale_color_gmri() +labs(y ="Temperature Anomaly\n(1991-2020 Baseline)",linetype ="",color ="Area")
Warning: Using the `size` aesthetic in this geom was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` in the `default_aes` field and elsewhere instead.
Warning: Removed 212 rows containing missing values or values outside the scale range
(`geom_point()`).
:::
Long-Term Trends and Net Change
These two datasets span a different range of years. For bottom temperature we will use from 1970 onward b/c we don’t have any other data before then in this project so we already are framing from 1970-onward anyways. Surface temperature will be done from 1982 onward because that is the first full year.
In the event of first order trends, net change over that time-span will be calculated as rate * years.
I’ll also be focusing on seasons to match the trawl survey as well.
Annual Trends
# Run the bottom tempoerature trends on annual summariesstart_yr <-1970; end_yr <-2019annual_bottom_trends <- annual_temps %>%filter(year %in%c(start_yr:end_yr)) %>%split(.$survey_area) %>%map_dfr(function(.x){ annual_mod <-lm(bot_temp ~ year, data = .x) mod_summ <-summary(annual_mod)[["coefficients"]]# summary info mod_trend =ifelse(mod_summ[2,4] <0.05, T, F) mod_rate =coef(annual_mod)[[2]]# tidy little output dataframe signif <-data.frame("season"=c("Annual"),"trend"=c(mod_trend),"rate"=c(mod_rate) )return(signif)},.id ="survey_area") %>%mutate(var ="Bottom Temperature\n(1970-2019)",temp_delta =ifelse(trend == T, rate * (end_yr - start_yr), NA))
# Run the bottom tempoerature trends on annual summariesstart_yr <-1982; end_yr <-2019annual_surface_trends <- annual_temps %>%filter(year %in%c(start_yr:end_yr)) %>%split(.$survey_area) %>%map_dfr(function(.x){# linear model annual_mod <-lm(sst ~ year, data = .x) annual_summ <-summary(annual_mod)[["coefficients"]]# summary info mod_trend =ifelse(annual_summ[2,4] <0.05, T, F) mod_rate =coef(annual_mod)[[2]]# tidy little output dataframe signif <-data.frame("season"=c("Annual"),"trend"=c(mod_trend),"rate"=c(mod_rate) )return(signif)},.id ="survey_area") %>%mutate(var ="Surface Temperature\n(1982-2019)",temp_delta =ifelse(trend == T, rate * (end_yr - start_yr), NA))
# Combine those jokersannual_trends_summary <-bind_rows(annual_surface_trends, annual_bottom_trends) %>%mutate(survey_area =fct_relevel(survey_area, rev(area_levels_long)))# Dumbell Plotggplot(annual_trends_summary) +geom_segment(aes(y = survey_area, yend = survey_area, x =0, xend = temp_delta, color = var), position =position_dodge()) +geom_point(aes(x = temp_delta, y = survey_area, color = var)) +geom_vline(xintercept =0) +scale_color_gmri() +scale_x_continuous(labels =label_number(suffix = deg_c)) +facet_wrap(~var, ncol =2) +theme(panel.grid.major.y =element_blank()) +labs(x =expression(Delta~"Temperature"), y ="", title ="Net Temperature Change",fill ="")
Warning: Width not defined
ℹ Set with `position_dodge(width = ...)`
# Maybe this should be a table:warming_table <- annual_trends_summary %>%select(-c(season, trend)) %>%rename(area = survey_area,`Net Change`= temp_delta, `Rate of Warming`= rate) %>%group_by(var) %>%gt(rowname_col ="area") %>% gt::tab_header(title ="Regional Surface and Bottom Temperature Change") %>%fmt(columns =`Rate of Warming`, fns =function(x){str_c(signif(x,2), deg_c, " / year")}) %>%fmt(columns =`Net Change`, fns =function(x){str_c(signif(x,2), deg_c)})warming_table
Are regional anomalies (either annual or monthly) too correlated with actual temperatures to include together in a model?
# Get overall correlationseasonal_temps %>%select(bot_temp, bt_anom) %>%drop_na() %>%summarise(corr =cor(x = bot_temp, y = bt_anom))
# A tibble: 1 × 1
corr
<dbl>
1 0.291
seasonal_temps %>%select(bot_temp, bt_anom) %>%drop_na() %>%ggplot(aes(bot_temp, bt_anom)) +geom_point() +labs(x ="Absolute Bottom Temperature", y ="Regional Bottom Temperature Anomaly", title ="Correlation: 0.29")
Modeling with Temperature & Anomalies
#### Load Data ####wigley_bmspectra_df <-read_csv(here::here("Data/model_ready/wigley_community_bmspectra_mod.csv"))
Rows: 398 Columns: 17
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (4): season, survey_area, spectra_type, area_titles
dbl (13): est_year, xmin_fit, xmax_fit, xmin_actual, xmax_actual, n, b, conf...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Joining with `by = join_by(survey_area, est_year, season)`
# What do they look like on the same x axis?# Add median weight here:wigley_medwt_df <-read_csv(here::here("Data/model_ready/wigley_community_medsize_mod.csv"))
Rows: 398 Columns: 16
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): survey_area, season, area_titles
dbl (13): est_year, n_species, numlen_adj_total, mean_len_cm, med_len_cm, mi...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
[conflicted] Will prefer lmerTest::lmer over any other package.
# Lets compare which one does betteractual_temps_model <-lmer( b ~ survey_area * season *scale(roll_temp) + survey_area*scale(log(total_weight_lb)) + (1|est_year), data = wtb_model_df)anomalies_temps_model <-lmer( b ~ survey_area * season *scale(roll_anom) + survey_area*scale(log(total_weight_lb)) + (1|est_year), data = wtb_model_df)# # Make it make sense# summary(two_temps_model)check_model(actual_temps_model)